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Abstract 

In this note, we introduce a new finite difference approximation called the Black- 
Box Logarithmic Expansion Numerical Derivative (BLEND) algorithm, which is 
based on a formal logarithmic expansion of the differentiation operator. BLEND 
capitalizes on parallelization and provides derivative approximations of arbitrarily 
precision, i.e., our analysis can be used to determine the number of terms in the 
series expansion to guarantee a specified number of decimal places of accuracy. Fur- 
thermore, in the vector setting, the complexity of the resulting directional derivative 
is independent of the dimension of the parameter. 

Key words: Finite difference algorithm, numerical differentiation, Taylor series 
expansions, sensitivity analysis, directional derivative 


1 Introduction 


Evaluation of derivatives is essential in optimization, control and sensitivity 
analysis. We consider the setting where the function of interest, 0:0—)- 
M, 0cR, is not available in closed-form but function evaluations are avail- 
able for any 6 G 0, i.e., we are in a “black box” scenario where derivatives 
of (j> have to be computed numerically by means of a finite difference (FD) 
approximation. For instance, the stationary distribution of a finite Markov 
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chain may fail to have a closed-form solution, which rules out analytical com- 
putation of derivatives, but the stationary distribution can be easily evaluated 
numerically. Many complex queueing networks fall into this category, and we 
use a simple tandem queue to illustrate our algorithm. 


The most well-known FD approximations, with h > 0, are the forward and 
central FD approximations given respectively by 


m 


<f>{0 + h) - (f)(9) 
h 


( 1 ) 


and 

_ <f>{0 + h) - <j>(6 - h) 
® ^ ~ 2 h 


( 2 ) 


Many other FD approximations involving weighted sums and differences of 
function evaluations on one-dimensional grids with arbitrary spacing are avail- 
able in the literature, including for higher-order derivatives; see, for example, 
[4], However, determining the weights required for each grid point is often 
computationally demanding and depends on the chosen grid. Also, most of 
the approximations involve recomputing the weights when additional preci- 
sion is required. 


In this note, we consider a simple FD series approximation that we call 
the Black-Box Logarithmic Expansion Numerical Derivative (BLEND), which 
is based on a formal logarithmic expansion of the differentiation operator. 
BLEND provides derivative approximations of arbitrary precision, in the sense 
that a specified number of decimal places of accuracy can be guaranteed by 
adding a sufficient number of terms in the series approximation, so the stop- 
ping criterion can be evaluated directly as the algorithm runs. The proposed 
BLEND algorithm capitalizes on modern computing platforms by exploiting 
parallel processing in evaluating the terms simultaneously. 


BLEND begins with the Taylor series expansion: 


oo in 

W + h) = £ -T VA>0. (3) 

,„o n - 

Denoting by A the differential operator A(f> := <\>' and by Th the shift-operator 
7h4> '■= 0(- + h), (3) can be rewritten in operator notation as 


W = £ 

n = 0 


h n 

n\ 


A> = £ 

n= 0 


(hA) n (j) 

n\ 


leading to the compact operator relationship Th = exp(hA). The existence of 
series expansions for inverses of analytic functions is a well-known result of 
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complex analysis (see [5,3]), and formal inversion yields 


A = 



( 4 ) 


for h sufficiently small, thus expressing the derivative in terms of the logarithm 
of the shift operator. 


Expanding the (natural) logarithm in (4) as a formal (Taylor) power series 
around the identity operator J gives the following formal representation: 


A 


U n = 1 n ' 

I T h ) n 

h Ti 

II n= 1 n 


( 5 ) 


for h sufficiently small. Noting that the shift operator Th satisfies 7)f = 

%h, k > 1 : 

(J- T h Y<t>=jz h)(-l)‘r^=E(-l + (6) 

for n > 1. Combining (5) with (6) yields the logarithmic series expansion of 
the derivative: 


W) = - 1 T £t— ("W + «0, (?) 

h n=\ n \ k J 

for h sufficiently small. A main contribution of this note is characterizing the 
domain of h for which (7) will return the correct answer, a result established 
in Section 3. 

This logarithmic series expansion idea leading to (7) was introduced by As- 
mussen and Glynn [1, pp. 212-213] as a potential basis for FD estimators in the 
Monte Carlo simulation setting. For example, using just the first term in the 
series (7) gives evaluation of two terms in (6) - specifically A « (Th — J)/h in 
shorthand operator notation of (5), leading to a “simulation-based first-order 
approximation” that is just the usual forward difference estimator correspond- 
ing to (1). However, [1, pp.213] concludes that due to numerical difficulties 
and variance challenges in the Monte Carlo simulation setting, “it is rare that 
estimators of higher order than the central difference estimator are used.” 

The main purpose of this note is to explore BLEND in the deterministic 
setting. Theoretical analysis and empirical evidence indicate that BLEND 
provides a practical alternative to other FD numerical approximations in many 
scenarios. In fact, the series expansion can also be used to derive unbiased 
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simulation-based finite- difference estimators, but these suffer from numerical 
instabilities similar to those observed in [1], 

The rest of this note is organized as follows. The main theoretical results for 
BLEND are given in Section 2, as well as the higher-dimensional extension 
to directional derivatives. The BLEND algorithm is introduced in Section 3. 
Numerical examples are presented in Section 4, and Section 5 concludes with 
a brief discussion of future application of BLEND to the stochastic setting. 


2 BLEND Analysis 


The main technical conditions required for BLEND are the following: 

(Cl) The function 0 is analytical on (6 — e, 6 + h 0 + e) for some e > 0. 
(C2) Fix N, and assume that for n < N: 

sup \^ n \§)\< Mb n , 

6<6<0+hn 

where 0^ = 0 and 0^ denotes the nth derivative of 0. 

The condition on the analyticity of 0 on (6 — e, 9 + ho + e) ensures that all 
higher-order derivatives are well defined on the interval of interest [6,6 + ho]. 
Note that if 0 is a polynomial of order k, then (C2) holds for b = max(h+hh, 1) 
and M — k\. Conditions (Cl) and (C2) allow a bounding of the n th power 
of the shift operator (J — 7h)" ■ 

Lemma 1 If (Cl) and (C2) hold, then for n < N : 

\(J -T,,) n m\< -f^{2hie) r ' . 

Vzt: n v 7 

Proof: For any n > N f and 0 < k < n, we have a Taylor’s series expansion 


0(h + kh ) — 0(h) + ^ (f>\6) 


( khf 
2 ! 


0"(h) + ...+ 


0 w (Sf) 


for some 6 {0,0 + kh). Multiplying each expression by ( — I j f, ( ? j and sum- 
rning up for k = 1 , ,n yields (below, ££ := 6 for 0 < i < n — 1 and any 
k < n): 
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(j - %rm = E(-i)‘u)^(« + kh ) 

=t(-D 

n ui n / \ 

=E|E(-1)‘*‘ RW) 

i=0 *' fc=0 W 

where, in the last equality, we have used the identity (Ruiz 1996, Corollary 2) 

E(- 1 )‘(^W = °. (8) 

valid for any polynomial V up to degree n — 1. Hence, for n > 1: 


\(j-T h rm\<M^±k^ n 

n! ^ 


k . 


„ r (2 nh) n n (2hbe) 
< M- — -— < M- J 


n\ \Z2nn ’ 

the last inequality following by Stirling’s approximation. 


(9) 


□ 


Example 1 Let f>(6,x) := dexp(— 0a;), a: > 0. Then for any n > 1 it holds 
that 


{J ~ %) n (j)( 6 , x ) 

= ^(-1)*Y ’’f] ( 9 + kh) exp [-(9 + kh)x] 
k= 0 \"V 

= dexp(— Ox) Yf-l) k ( ^ ] exp(— khx) 
k = o \"v 

+Lexp(— ^(— ?? J exp(— khx) 

k i v A v 

= dexp(— 0x) [1 — exp(— for)]” 

—nhex p [— (0 + /i)x] [1 — exp(— /uc)] n_1 . 

Dividing by —hn and summing everything up for n > 1, yields 

_jY — — — = (1 - frr) exp(-frr) = ^(0,x), (10) 

h n > 1 n <w 

which proves the validity of (5) for 4>(-,x). 
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Regarding the conditions of Lemma 1, we note that if 0 is an analytic function 
defined only on a bounded interval (/, b), then one can use a change of variables 
u : (a, oo) — * (/, b ) to obtain a new function p = 0 o u : (a, oo) — » M satisfying 
condition (Cl) and then recover the derivative 0'(0) by evaluating <// at the 
point u^iQ). Furthermore, recall that in general, analyticity amounts to the 
fact that for any compact set D, there exists M, b > 0 such that 


sup 

eeD 


(n) 


(*) 


< Mb n n\. 


In that sense, condition (C2) is more restrictive, as it imposes certain growth 
rates on increasing intervals. Under the conditions of Lemma 1, we see that 

h n 

0(0 + h) = y 0(»)(0) < Mexptbh), 

Z>, d n! 


which suggests that the result applies to exponentially bounded analytic func- 
tions. 


We conclude this section by discussing the extension to directional derivatives. 
Let 0 e and assume that 0 is analytical as a mapping of 0. The directional 
derivative of 0 is direction v := (n 1; . . . , v m ) is given by 


V*0(0) 


0(0 + hv) — 0(0) 

Inn — — 

h\\v\\ 


E 





where ||n|| denotes the Euclidean norm. Evaluating a directional derivative 
via the above analytical approach requires the evaluation of m partial deriva- 
tives. Alternatively, we can apply the shift-operator and obtain the following 
logarithmic series expansion of the directional derivative 

+ f=v). (ii) 

Note that the complexity of the logarithmic series expansion is independent 
of the dimension m of the parameter. 


3 The BLEND Algorithm 


We consider the finite truncation of the sum in (5) as an approximate value 
for the derivative, i.e., the BLEND approximation of order N and difference 
h for the derivative of 0(0) with respect to 0 is given by 

1 N n (—l) k f n \ 

a(m, h)m = -v e E 00 L#+ kh < 12 ) 

n n=lk=0 U \ K ) 
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and the error of the BLEND approximation will be denoted by 


R(N,h,m) = \^m - 


In the following we show how Lemma 1 can be used to determine the trunca- 
tion index N such that at least the first n digits of the derivative are exact. 


Lemma 2 Under conditions (Cl) and (C2), 


R{N,h,<t>{6 )) < 


M (2 hbe)( N+ V 

V2tt(N + 1)^ 1-2 hbe ’ 


N>1, 


provided h < 1/2 be. 


Proof: It follows from Lemma 1 that 


Vn> 1: \(J-T h ) n m\<M 


(2 hb 2 e) n 
y/2wn 


Hence, for h > 0 


E 

n=N + 1 


(j - r h ; 


n 


< 


M 


(N + 1) 2 \/2n n= N+i 
M 


E (2 hbey 


(N+l)*y/2w 


(2 hbe] 


N+l 


1 — 2 hbe ’ 


provided h < l/2be. 


□ 


Example 2 Consider the mapping (f)(9) = sin(0). We apply BLEND for 
computing the derivative of sin(d) at 9 = 0. Condition (C2) holds for all 
N with M = 1 = b. By Lemma 2, BLEND yields the correct result for 
h < l/2e ~ 0.1839. As can be seen in Table 1 executing the BLEND algo- 
rithm for N = 8 yields 8 exact digits of the derivative. To illustrate the impact 
on the bound on h put forward in Lemma 2, we present in Table 2 the output 
of BLEND fork = 1. 


As can be seen in Table 2 executing the BLEND algorithm for h — 1 fails to 
produce the correct output, which stems from the fact that 1 > l/2e ~ 0.1839. 

The fact that BLEND fails to yield the correct output for large h can best be 
seen by applying BLEND to the derivative of sin(d) at 9 = 0 with h = 2n, 
which yields 0 = A (f>(N, 2i r) for all N, whereas 1 is the correct answer. 

Lemma 2 shows that the error of considering A (N, h)<f>(9) rather than A (h)<f>(9) 
is of order 0((2/?he) ( ^ Ar+1 ^). Moreover, for given N, provided h is sufficiently 
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N A (IV, 0.01) 

1 0.998334166468282 

2 1.003321678961257 

3 1.000029893016725 

4 0.999980308400858 

5 0.999999646316608 

6 1.000000137620388 

7 1.000000003815154 

8 0.999999998963623 

true 1.0 

Table 1 

BLEND Approximation for the derivative sin(0) at 6 = 0, small h 

N A(1V, 1.0) 

1 0.841470984807897 

2 1.228293256202952 

3 1.207506816871789 

4 1.015352293328013 

5 0.885486080979581 

6 0.903764738896000 

7 1.003453862663737 

8 1.071046882890327 


Table 2 

BLEND Approximation for the derivative sin(0) at 9 = 0, large h 

small so that it satisfies the condition in Lemma 2 and upper bounds for M 
and b are known, solving for h in 


M (2 hbe) N+1 
1-2 hbe 


10 -(K+1) 


(13) 


yields a value for h such that the FWD (forward finite difference) approxi- 
mation A (N, h)4>{9 ) is exact in at least the first K digits. We call this the 
K -exact FWD approximation. 


Typically, M and b cannot be computed exactly. In this case the bound on 
the remainder put forward in Lemma 2 can be facilitated through the ge- 



ometric convergence rate of the series in (12) for small h. Specifically, in- 
specting the values of A (IV, h ) for h fixed and N = 1, 2, . . . one expects that 
after a possible transient phase the geometric error decrease established in 
Lemma 2 sets in. This gives rise to the following heuristic: Compute A (IV, h), 
for N — 1,2,..., fV max , with iV max a pre-specihed number. If no stabilization 
is seen, then decrease h and compute again the A max values. If stabilized with 
the first L digits being the same for 7V max _ j and fV max , then in light of the 
geometric error decrease, accept the first L digits in A (IV', h) - as if the values 
of A (IV, h) for N > N' only affected digits L + 1 and larger. We call this the 
stabilization indication of the series. We summarize this in BLEND algorithm 
presented in the following. 


BLEND Algorithm 

(%) Choose h small and set 7V max . 

(ii) Evaluate A (IV, h)(j>(9) in (12) for N — 1, , fV max . 

(in) Let L be such that the first L digits in A(lV max — 1, h)<f>(0) and A(fV max , h)cj)(6) 
are equal, then the value of A(fV max , h)4>(9) up to first L digits is outputted 
as the exact value for the first L digits of the derivative. 

Remark 1 The BLEND exploits parallel processing by executing each of the 
evaluations in step (ii) simultaneously. As an illustration, 7V max = 8 would 
be a natural choice for many current portable multi-core computing platforms 
that are commonly available, e.g., Apple laptops. 


4 Numerical Examples 


We illustrate the K - -exact FD approximation as well as the BLEND algorithm 
for a series of examples, beginning with some toy examples and then consid- 
ering a more complicated example that illustrates the type of setting that we 
envision for the main application of BLEND. For comparisons, recall that the 
FD approximation is just the first term in the tables below, when N — 1. 

Example 3 Let f>(9) = 5 9 4 and note that (Cl) holds for ho = oo. We will 
compute the derivative at 9 0 = 2. In order to apply condition (C2) , we assume 
that h is no larger than h = 0.1. Only the first f derivatives are significant 
and we arrive at 


sup 

0o<0<0o-\-hn 


(n)l 


(9) < 120(2 + 0.1 x 4) 4 , 1 < n < 4, 


and thus M = 120 and b = 2.4. 
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We apply BLEND for N max = 8 and, h = 0.01. The numerical results are 
provided in Table 3. The numerical results put forward in Table 3 show that 


N 

A(N, 0.01) 

1 

160.1200400049834 

2 

159.9999199699909 

3 

160.0000299999870 

4 

159.9999999999799 

5 

159.9999999999719 

6 

159.9999999999577 

7 

159.9999999999281 

8 

159.9999999999981 

true 

160 


Table 3 

The Finite Approximation for the derivative 5 0 4 at 9 = 2, h = 0.01 

our stabilization indication works in natural way. Notice that up to 10 digits 
A(8,0.01) is rounded off to 160, giving same digits up to precision 10~ 10 . 

We now turn to K -exact FWD approximation. Suppose we want to compute 
the first 6 digits of the derivative exactly, with lV max = 2. Then we solve (13) 
for h with M = 120 and b = 2.4, i.e., 

120 (4.8 hef = 7 

v /2^25 1 - 4.8 he 

so that h has to be smaller than 0.00013. The resulting value for A(2, 0.00013) 
is 1.599999999999915. 

In the following we illustrate the application of our results to directional deriva- 
tives. 

Example 4 Consider 0 € M m and let 


4 ( 0 ) = 

i= 1 

for some constants a* G M, 1 < i < m. We consider the directional derivative 
of (f)(9) in direction v, where v e [—1, l] m is chosen such that Jf™=i v i — 0. 

For our numerical experiment, we let m = 9, a* = 2~\ for 1 < i < m, and 

t7=(l/3)(-l, 1,-1, -1,1, 1,1, -1,1), 
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so that ||u|| = 1. We obtain 

d( t ) /m O t 

= for\<i<m. 

We now apply BLEND for iV max = 8 and h = 0.01 for computing the direc- 
tional derivative of (f)(9) in direction v at 9, with 6i = i, for 1 < i < m. The 
numerical results are provided in Table f. 


N 

A(JV, 0.001) 

1 

3.958029296875054 

2 

3.957031250000576 

3 

3.957031250002056 

4 

3.957031250004276 

5 

3.957031250005520 

6 

3.957031250007444 

7 

3.957031250013154 

8 

3.957031250013043 

true 

3.9570312500138101 


Table 4 

The BLEND Finite Approximation for the directional derivative 

Now we consider the queueing example, where the constants required to set 
the K -exact FD approximation would not be available. 

Example 5 Consider a two-station tandem queueing system with finite ca- 
pacity Ni at station 1 and N 2 at station 2. Jobs arrive to the network accord- 
ing to a Poisson process with arrival rate X, and service times are independent 
and identically distributed (i.i.d.) exponential with rate /i, at station i — 1,2. 
When there is no waiting place available at station 1, an arrival is rejected 
and lost. When there is no waiting place available at station 2, service station 
1 is stopped. This example is taken from [2], and we refer for motivation and 
details to the references therein. Due to the finite buffers, no closed-form ex- 
pression for the stationary distribution exists. However, letting Q denote the 
infinitesimal generator of the process, the stationary distribution solves nQ = 0 
with normalizing equation — 1 and is easily numerically evaluated. 

We apply the BLEND algorithm for computing the derivative of the blocking 
probability with respect to X, taking X — 1, fi — 1, r/i = 1 and r/ 2 = 2, with 
finite capacity queue sizes N { — N 2 — 10. The numerical results are given in 
Table 5, where the true value has been obtained by the finite difference method 
through a series of experiments. 
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N A (IV, 0.01) 


1 

0.613180514116096 

2 

0.610046682208255 

3 

0.609671969013386 

4 

0.609661671019043 

5 

0.609662935724646 

6 

0.609663162694883 

7 

0.609663173459084 

8 

0.609663170509458 

true 

0.609663168 


Table 5 

The BLEND Approximation for the Loss Probability Sensitivity with respect to A 
at A = 1 

5 Conclusion and Future Research 


We presented a new finite difference derivative approximation called the BLEND 
algorithm. BLEND is particularly useful when the expression of 0 (the func- 
tion of interest) is not available in closed-form, but the values at arbitrary 
points within its domain can be efficiently numerically computed. We charac- 
terized the value of the difference parameter h for which the BLEND algorithm 
applies, and provided a bound on the error that leads to an estimate on the 
number of terms required to achieve a particular degree of precision. Various 
numerical examples illustrate the effectiveness of the BLEND algorithm, in- 
cluding the importance of choosing the parameters correctly and the practical 
implementation using parallel computing. 

Future research of interest is considering the extension to the stochastic case 
by letting 0(d) = E [Z(9)\, for Z(9) the underlying stochastic variable, and the 
application of BLEND to high dimensional problems. 
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